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Abstract 

We perform a complete theoretical analysis of the azimuthal angular correlation of two-hadron 
productions in the forward dAu collisions at RHIC in the saturation formalism, and obtain a 
very good agreement with the experimental data. It is demonstrated that the suppression and 
broadening of the away side peak provide a unique signal for the onset of the saturation mechanism 
at small- cc in a large nucleus. We emphasize that future experiments of di-hadron correlations in 
pA collisions at both RHIC and LHC, and in eA collisions at the planned electron-ion collider, 
shall provide us with a thorough study and understanding of the strong interaction dynamics in 
the saturation regime. 
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Introduction. Two particle correlations in high energy hadronic scattering processes have 
played important roles in studying the strong interaction QCD dynamics jlj, [2j, novel phe- 



nomena in proton-proton collisions 



6]. 



3| , and the medium effects in heavy ion collisions |4 
In particular, the azimuthal correlation of two particles with large transverse momenta pro- 
vide critical evidence of the strong jet quenching effects in the hot dense medium created in 
the heavy ion collisions at RHIC and the LHC. An important advantage of these measure- 
ments is that they reveal the physics by the observables themselves, which do not need, e.g., 
the proton reference to quantify the phenomena. It was suggested in Ref. to study the 
cold nuclear matter in dense region by measuring the forward di-jet (di-hadron) production 
in proton-nucleus collisions, where the de-correlation of the back-to-back dijet in the forward 
pA collisions can be used to signal the gluon saturation at small-x in a large nucleus. This 



prediction was qualitatively confirmed by the STAR 



[9] 



and PHENIX M collaborations from 



di-hadron correlation measurements in dAu collisions at RHIC, and have been considered 
as the best evidence for saturation physics. Early attempts [h], have been made to 
understand the experimental data quantitatively, where, however, inappropriate approxi- 
mations were taken in the calculations. In particular, in this paper we correctly distinguish 
the Weizsacker- Williams (WW) gluon distribution from the dipole gluon distribution. In 
this work, we will provide, for the first time, a quantitative and thorough description of the 
experimental data in the saturation formalism, including the large broadening of the angular 
distribution and suppression of the peak for the away-side hadron. Our results emphasize 
the relevance of the gluon saturation in the kinematic regions covered by the STAR and 
PHENIX collaborations and signify the usage of the two-hadron correlation as an important 
tool to investigate the QCD dynamics in the small- x limit. 



There have been theoretical arguments 12l4l5{| which suggest that the gluon distribution 



saturates at small Bjorken-x. The color-glass-condensate model (CGC), has been proposed 
to describe the gluon saturation phenomenon at small-x. An important feature of this 
approach is the appearance of the dynamically generated saturation scale Q s (x) which sep- 
arates the dilute and dense partonic regimes. The experimental data from HERA are well 
described with the saturation model calculations, with the saturation scale of the order of 
~ 1 GeV. Furthermore, the forward hadron production in nucleon-nucleus (pA) collisions 



has been systematically studied 16) in the CGC formalism, where the unintegrated gluon 



distributions (UGDs) are important ingredients to describe the phenomena. They unveil the 



importance of the multiple interaction effects in the factorization of the hard processes in 
the small- x calculations. Nevertheless, there exists much more interesting dynamics [l^] in 
saturation physics which can only be explored by di-jet or di-hadron production processes 
as we will demonstrate in the following calculations. 

In this paper, we focus on two-particle production in the forward direction of pA (dAu 
at RHIC) collisions, 

p + A^h x + h 2 + X , (1) 

where two hadrons hi and hi with large momenta are produced. The above process is 
sensitive to the gluon distributions at small- x in the nuclear target. In order to correctly take 
into account the multiple interaction effects, we follow the CGC framework to calculate the 

n 

two particle production [17[ . An effective k t factorization can be established for this process 
in the back-to-back correlation limit, and the differential cross sections can be expressed 
in terms of various UGDs, which can be related to two fundamental UGDs: the dipole 
gluon distribution xG^ 2 \x, q±), and the WW gluon distribution xG^(x,q±). Only with 



this effective kt factorization, can one describe a. 



and suppression) of the STAR 



and PHENIX 



lthe features (including both broadening 
data systematically. These results also 



agree with previous calculations for two-particle production in pA collisions in the general 



kinematics region 18|, [19] 



In the RHIC experiments, the di-hadron correlations are measured by the coincidence 
probability C(A<f>) = A pail .(A0)/A^ rig , where iV pair (A0) is the yield of two forward tt° which 
includes a trigger particle with a transverse momentum p^[ s and an associate particle with 
pff an d the azimuthal angle between them A0. We calculate the single and two-particle 
cross sections and obtain, 



Ji 



dcrP A ^ h ± h 2 



J\pi±\ dy 1 d' 2 p 1± 

where the dependence on the rapidities of the two particles is implicit. 

Single inclusive cross section. Let us first discuss the single inclusive hadron production. 



The leading-order single inclusive cross section 16) in pA collisions is given by the product 



of the integrated parton distributions of the projectile proton and the unintegrated gluon 



distributions of the target nucleus: 

d 2 bd 2 P± dy h = J Zh H 

+ x p g f (x p )F Xg (k ± )D h/g (z 1 ) , (3) 

where the sum over quark flavor is implicit, b represents the impact parameter in pA colli- 
sions, p± and yh are transverse momentum and rapidity of the hadron, q(x p ) and g(x p ) are 
integrated quark and gluon distributions from the projectiles, D(z) the associated frag- 
mentation functions with p± = Z\k±, x p = p±e m / 'z\^/s and x g = p±e~ Vh / Ziy/s. The 
dipole gluon distributions F x (k±) and F Xg (k±) are Fourier transform of the dipole scatter- 
ing amplitude in the fundamental and adjoint representations, respectively. In particular, 
F Xg (k±) oc x g G^ 2 \x g , k±)/k 2 L . In terms of the numerical study, we are able to describe the 
forward single hadron production cross sections measured by both BRAHMS and STAR 
up to p± = 3.0 GeV with a K-factor about 0.8 for y h = 2.0 and 0.5 for Uh = 3.2. In this 
numerical evaluation, we follow the NLO sets of MSTW parametrizations [20| for the parton 
distributions and DSS parametrizations 21 1 for the fragmentation functions 

Two-particle production in forward pA collisions. Two-particle production contains the 
correlated and uncorrelated contributions, 

da { P A->h lh2 ) = da ( p A-* hl h 2 ) + _ ( 4 ) 

The correlated hadron production comes from the partonic 2 — » 2 processes, where these two 
particles are back-to-back correlated and form the away side peak in the azimuthal angular 
distribution (A0 = 7r). The near side correlation comes from the particle decay or the 
same jet fragmentation if they are at the same rapidity. In this letter, we will focus on the 
back-to-back correlation region, namely the away side peaks. According to Ref. [171 ]. we can 
write down the differential cross section for the two-particle production in the back-to-back 
correlation limit, 

daiorr7 fhlfl2 ^ f dzi dz 2 a 2 s 
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y 2 ~2 £2 



[x p q(x p ) J-qg 



dy hl dy h2 d 2 p 1± d 2 p 2 ± J z\ z\ 
x H gg {D hl/q ( Zl )D h2/g {z 2 ) + D h2/q ( Zl )D hl/g (z 2 )) 

+x p g(x p )F®H®D hl/g (z 1 )D h2/g (z 2 )] , (5) 



1 A recent next-to-leading order calculation for inclusive hadron production suggests that the appropriate 



choice for the factorization scale to be around the saturation scale 
our calculations. 



22[. We have followed this choice in 



4 



where x g = X\e V1 + x 2 e m and x p = x\e Vl + X2e y2 with X{ = \ki±\/y§ and k{± = Pi±/zi, 
and if « are various UGDs and the associated hard coefficients, respectively. Their 
expressions can be found in Ref. 17| . The partonic center of mass energy squared s is 
defined as s = P 2 J z(l - z) with P ± = (k x± - k 2± ) /2 and z = x l e^/x 2 e y \ In the CGC 
calculations [17j|, P± = (1 — z) ku_ — zk 2 ± also enters in the hard coefficients, which equals 
to P± in the correlation limit. The difference between P± and P± will be used to estimate 
the theoretical uncertainties in the following calculations. In the typical kinematics of the 
forward collisions at RHIC, we find that x p ~ 0.1 and x g < 10~ 3 , where both the quark 
initiated processes (q — > qg channel) and gluon initiated processes (g — > gg) contribute. 



Comparing the above equation to Eq. (3) of Ref. 10J, one immediately finds notable 



differences between the results. In particular in Ref. [10[ the only channel calculated was 



q — > qg. Moreover, in this channel our results do not agree with results in Ref. [10| since 
in the latter work the contributions from the WW gluon distribution were not taken into 
account. These contributions are essential in order to reproduce correctly the collinear 
factorization results for dijet production in the dilute limit. 

The unintegrated gluon distributions in Eq. fl5]) are largely un-explored, in particular, 
for those related to the WW gluon distribution. The energy evolution is important to un- 
derstand their behavior depending on x g , of which for the dip ole gluon distribution, the 



Balitsky-Kovchegov (BK) evolution, has been well studied 15( and demonstrated the so- 



called geometric scaling 123] in the solution. The scaling was found to be related to the trav- 



eling wave solutions [24 125| of the BK evolution. The energy evolution equation for the WW 
gluon distribution has recently been systematically investigated [26[]. An important result 
from these studies is the geometric scaling similar to the dipole gluon distribution. There- 
fore, as a first step, we can parametrize these gluon distributions from a model calculation, 
and include the energy dependence by assuming the geometric scaling and x s -dependence of 
the saturation scale. In the following, we adopt the Golec-Biernat Wusthoff model [^j for 
the dipole gluon distribution which successfully describes the low-x DIS structure functions 



at HERA, then extend it to the WW 
by modifying the saturation scale as 



gluon distribution and include the nuclear dependence 



15| 



Q'a = c{b)A^Ql( x ) , (6) 
where c(b) represents the profile function of nucleus depending on the impact parameter b of 
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ixlx )- x with Q s0 = 1 GeV, x = 3.04 x 1(T 3 and A = 0.288 follow 



the collision, Q 2 s {x) = Q 2 sQ 
GBW parameterizations 27|. The profile function c(b) is closely related to the centrality 
of the pA (or dA) collisions. Central collisions give large value of c(b), while peripheral 
collisions correspond to small profile function. 

We would like to emphasize that the GBW model is not sufficient to describe the UGDs 
in the region that k± is much larger than Q s . However, for the forward pA collisions, the 
saturation scale Q s is large enough to cover most of the kinematics where k± is around 
Q s and we will be able to well describe the experimental data. For pp collisions, we have 
to either modify the GBW model or include the broadening effect from the fragmentation 
function to describe the experimental data. 

Double parton scattering contribution. Now, we turn to the un-correlated two-particle 
production in the process ([T|). This part mainly comes from two independent hard scat- 
terings, which is referred as double parton scattering (DPS) contributions (see recent de- 



velopments |28l432j]). It has been pointed out in Ref. |28[ that the DPS may exceed the 
single parton scattering contribution in the forward pA collisions. Following these ideas, 
we estimate its contributions in pA collisions in the saturation formalism. In particular, 
the multiple interactions from the nuclei side has been taken into account in the CGC fac- 
torization formalism [l^. For the proton side, we follow a simple parametrization for the 
double parton distribution: V l J(x p ,x' p ) = C(x p ,x' p )fi(x p ) x fj(x' p ) with C ~ 1, where i and 
j represent the two partons from the nucleon which participate the hard scattering, x p and 
x' p for their momentum fractions. The final expression reads as 

^CTuncorr. 1 ^ f dZ\ dZ2 



* £x p /^K/ 3 (4)Fi-HM^f(fc2±) ■ (7) 

ij 

where ij represent flavors of the partons from the nucleon and the associated UGDs from 
nuclei, and x p and x g are determined by the kinematics of the two hard scatterings. 

An important feature of the above DPS contribution is that the two hard scatterings are 
independent to each other at the leading order approximation 



28 



301 ] . Therefore, the two 



particles in the final state are un-correlated, and their azimuthal angle distribution will be 
flat. This leads to the so-called pedestal contribution in the experimental measurements. 

Comparison with the experimental data. With the above formulas, we are ready to com- 
pare to the experimental data on the two-particle correlation measurements in the forward 
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FIG. 1. The forward di-pion correlations C(A<p) of Eq. ([2]) at y\ ~ yi ~ 3.2 in peripheral and central 
dAu collisions compared to the preliminary data from the STAR collaboration [8]. Centrality 



definition follows Ref. 



8J], where the average impact parameters are found around 6.7fm and 2.7fm 



accordingly. The grey error band comes from using P± or P± in the hard coefficients in Eq. ([5]). 



dAu collisions at RHIC. Before we do that, we would like to emphasize that the saturation 
scale plays a key role in describing the correlation C(A0) of the away side peak, including 
both broadening and the suppression. First, the width of the away side peak will increase 
with the saturation scale because of the broadening effects. Quantitatively, the effective 
fc t -factorization formula of Eq. fl5]) lead to stronger broadening effects compared to the naive 
fc t -factorization calculations. This is because the various gluon distributions contain the 
convolution of the UGDs and will enhance the broadening. Without this enhancement, we 
can not describe the broadening effects. In particular, when the saturation scale reaches the 
transverse momenta of the dijet, the away side peak will almost disappear as indicated in 
the experimental data for the central collisions at RHIC and the theory calculations as well. 

Second, the magnitude of the correlation C(A<f>) is also sensitive to the saturation scale 
Q s . In particular, larger Q s push the dipole gluon distribution to larger transverse mo- 
mentum, which leads to single particle production fl3]) increasing with Q s . The correlated 
two-particle production cross section (JSJ), however, decreases with Q s for the same reason. 
Therefore, the correlated contribution to C(A0) decreases accordingly. Our numeric eval- 
uation also supports this conclusion. On the other hand, the un-correlated two particle 
production cross section ((7j) roughly depends on the product of two single particle cross sec- 
tions. Therefore, its contribution increases more rapidly with Q s than that of single particle 
cross section. The consequence is that the pedestal contribution increases with Q s . 



All these features are evident when we compare to the STAR data 8]. As an example, we 
show in Fig. 1 the results for p*" g > 2 GeV and 1 GeV < p^f < p^ at y\ ~ 1/2 ~ 3.2 in the 
peripheral and central collisions, respectively. In our calculations, we assume a fixed strong 
coupling constant a s = 0.35. The saturation scale Q s (c(b) in Eq. ([6])) is the only parameter 
to fit the data, for which we found c(b) = 0.45, 0.56, 0.85 for the peripheral, minimum bias, 
and central collisions, respectively. These parameters are consistent with c(0) ~ 0.9 and 
the centrality dependence of the nuclei profile for these collisions, by using either the hard 
sphere model c(b) = c(0)a/1 — b 2 / R 2 A or the Wood-Saxon model. The pedestal contributions 
to C(A0) are found to be around 0.016 and 0.018 for the peripheral and central collisions, 
which are in rough agreement with the experimental measurements. In the plots of Fig. 1, 
in order to better compare the results, we used the experimental extractions of pedestal 
contributions. The grey error band of the theory calculations comes from the difference 
between P± and used in the hard coefficients in Eq. fl5]). Similar results are also found 
when we compare to the correlation measurements from PHENIX collaboration. 

It is important to note that in the central dAu collisions the disappearing of the away-side 
peak in the di-hadron production indicates that the saturation scale is the same order as 
the hard probe of the jet transverse momentum, which is a clear signal of the onset of the 
saturation mechanism for this observable. From the kinematics of this collision, we conclude 
that the saturation scale Q s reaches at ~ 2 GeV at x g ~ 6 x 10~ 4 in the center of the 
gold nucleus with jet transverse momentum k±_ ~ 3 GeV at rapidity 3.2. We hope that the 
future measurements at RHIC and LHC will provide more information and help to map 
out the complete phase structure of the cold nuclei at small- x . We want also to emphasize 
that the above conclusion is very general and independent of the model we used for the 
UGDs. As mentioned above, the GBW model captures the main features for the UGDs at 
the transverse momentum around the saturation scale in which most of the data exist. 

Meanwhile, the PHENIX collaboration has also reported the nuclear suppression factor 
JdA defined as 

1 pair / 

j _ 1 <? d A l°dA / q\ 

dA ~ /AT \ Pair / ' ^ 

\iV coU ; ap p jOpp 

where cr pair is the cross section of two-particle production in dA and pp collisions with the 
pedestal contributions subtracted. a d A and the cross sections for the full event 

selection. In the absence of gluon saturation and nuclear effects, the dihadron cross-sections 
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FIG. 2. The nuclear suppression factor J^a of Eq. ([8j) for the two particle production in dAu 
collisions as function of Xj^f. The experimental data from PHENIX [9] and theory calculations 
are pedestal contributions subtracted. 



are expected to scale with (iV coll ). Therefore, J^a should be equal to unity in the dilute 



regime and suppressed in the dense regime. The experimental kinematic variable 



(pi±e~ yi + P2J_e~ 2/2 )/-y/i has been used to represent the ^-dependence [9]. In Fig. 2, we 
calculated J ha as function of xj^ for a typical transverse momentum pu_ = p2j_ = 1.0 GeV 
with the same rapidity for the two particles, for central collision (c(b) = 0.85) and peripheral 
collision (c(b) = 0.45), where we have used P± in the hard coefficients in Eq. (jHJ). Similar 
results are obtained with choice of P±. As a comparison, we list the PHENIX data in Fig. 2, 
where the data are for different values of Pi± and rapidities. From this figure, we clearly see 
that the the suppression of J^a at low a; A * s due to the saturation of the cross sections in 
pp and gL4 collisions, as also indicated by the experimental data. However, we emphasize 
that JdA depends on the pp reference, for which our model calculations should be taken with 
cautions. 

In summary, we have carried out a complete numerical study on the forward dihadron 
correlations in dAu collisions and found good agreement with the experimental data from 
RHIC. These results demonstrated that the saturation formalism developed recently can be 
used to describe the broadening and suppression of the away side peak in the di-hadron 
production in pA collisions. This emphasizes that the di-hadron (dijet) correlation provides 
a unique signal for the onset of saturation mechanism at small- x in a large nucleus. Future 



experiments at RHIC, EIC and LHC will provide excellent opportunities to thoroughly 
investigate the QCD dynamics in the saturation regime. 

We thank E. Avsar, L. Bland, M. Chiu, F. Dominguez, C. Marquet, L. McLerran, A. 
H. Mueller, J.W. Qiu, M. Strikman, R. Venugopalan, W. Vogelsang and N. Xu for helpful 
conversations. In particular, we are grateful to R. Venugopalan for stimulating discussions 
and comments. This work was supported in part by the U.S. Department of Energy under the 
contracts DE-AC02-05CH11231 and DOE OJI grant No. DE - SC0002145. We are grateful 
to RIKEN, Brookhaven National Laboratory and the U.S. Department of Energy (contract 
number DE-AC02-98CH 10886) for providing the facilities essential for the completion of this 
work. 



[1] V. M. Abazov et al. [ DO Collaboration ], Phys. Rev. Lett. 94, 221801 (2005). 

[2] V. Khachatryan et al. [ CMS Collaboration ], Phys. Rev. Lett. 106, 122003 (2011). 

[3] V. Khachatryan et al. [ CMS Collaboration ], JHEP 1009, 091 (2010). 

[4] J. Adams et al. [ STAR Collaboration ], Phys. Rev. Lett. 91, 072304 (2003). 

[5] G. Aad et al. [ Atlas Collaboration ], Phys. Rev. Lett. 105, 252303 (2010). 

[6] S. Chatrchyan et al. [ CMS Collaboration ], Phys. Rev. C84, 024906 (2011). 

[7] C. Marquet, Nucl. Phys. A 796, 41 (2007). 

[8] E. Braidot, for the STAR Collaboration, Nucl. Phys. A 854, 168 (2011). 

[9] A. Adare et al. [PHENIX Collaboration], Phys. Rev. Lett. 107, 172301 (2011). 
[10] J. L. Albacete and C. Marquet, Phys. Rev. Lett. 105, 162301 (2010). 
[11] K. Tuchin, Nucl. Phys. A846, 83-94 (2010). 

[12] L. V. Gribov, E. M. Levin and M. G. Ryskin, Phys. Rept. 100, 1 (1983). 
[13] A. H. Mueller and J. w. Qiu, Nucl. Phys. B 268, 427 (1986). 

[14] L. D. McLerran and R. Venugopalan, Phys. Rev. D 49, 2233 (1994); Phys. Rev. D 49, 3352 
(1994). 

[15] F. Gelis, et al, Ann. Rev. Nucl. Part. Sci. 60, 463-489 (2010); and references therein. 
[16] A. Dumitru and J. Jalilian-Marian, Phys. Rev. Lett. 89, 022301 (2002); J. L. Albacete, C. Mar- 
quet, Phys. Lett. B687, 174-179 (2010). 
[17] F. Dominguez, B. W. Xiao and F. Yuan, Phys. Rev. Lett. 106, 022301 (2011); F. Dominguez, 

10 



et al, Phys. Rev. D 83, 105005 (2011). 
[18] J. P. Blaizot, F. Gelis and R. Venugopalan, Nucl. Phys. A 743, 57 (2004). 
[19] J. Jalilian-Marian and Y. V. Kovchegov, Phys. Rev. D 70, 114017 (2004) [Erratum-ibid. D 

71, 079901 (2005)]. 

[20] A. D. Martin, W. J. Stirling, R. S. Thorne and G. Watt, Eur. Phys. J. C 63, 189 (2009). 
[21] D. de Florian, R. Sassot and M. Stratmann, Phys. Rev. D 75, 114010 (2007); Phys. Rev. D 

76, 074033 (2007). 
[22] G. A. Chirilli, B. -W. Xiao and F. Yuan, arXiv:1112.i06T1 [hep-ph]. 
[23] A. M. Stasto, K. J. Golec-Biernat and J. Kwiecinski, Phys. Rev. Lett. 86, 596 (2001). 
[24] S. Munier and R. B. Peschanski, Phys. Rev. Lett. 91, 232001 (2003); Phys. Rev. D 69, 034008 

(2004). 

[25] A. H. Mueller, D. N. Triantafyllopoulos, Nucl. Phys. B640, 331-350 (2002). |h ep-ph/0205167] . 
[26] F. Dominguez, et al, Phys. Lett. B 705, 106 (2011); A. Dumitru, et al, Phys. Lett. B 706, 

219 (2011); E. Iancu and D. N. Triantafyllopoulos, JHEP 1111, 105 (2011). 
[27] K. J. Golec-Biernat and M. Wusthoff, Phys. Rev. D 59, 014017 (1998). 
[28] M. Strikman and W. Vogelsang, Phys. Rev. D 83, 034029 (2011). 
[29] J. R. Gaunt, et al, Eur. Phys. J. C 69, 53 (2010); and references therein. 
[30] M. Diehl and A. Schafer, Phys. Lett. B 698, 389 (2011). 
[31] J. R. Gaunt, W. J. Stirling, JHEP 1003, 005 (2010). 

[32] B. Blok, Y. Dokshitzer, L. Frankfurt, M. Strikman, Phys. Rev. D83, 071501 (2011). 



11 



